%%%Main File for BCDR Simulations - gains from industrial policy and trade
%%%policy

%This file computes the gains from a given trade and industrial policy
%using the SOE assumption

%Version with intermediate goods and capital, economies of scale in
%gross inputs

%CES Version

clear all;

K=32; %number of sectors
N=62; %Includes ROW
K_man=15;    %number of manufacturing sectors

%%%Load data
%%%Data is organized as follows: first column origin country codes, second column
%%%sector codes, third column sectoral value added, fourth column aggregate origin trade balance,
%%%fifth column origin-sector expenditurs, 6th column starts the bilateral
%%%trade flows

A=importdata('simulation_ge_data_cap_int_avg_io.csv');
B=importdata('sim_params_cap_int_go.csv'); %these are the estimated elasticities

%Data vectors
X_ink = A(:,7:6+N); %Matrix of trade flows, (NxK)xN
D_i=A(1:N,4); %vector of trade balances, Nx1
GO_ik=A(:,3); %vector of sectoral gross output, (NxK)x1
expenditures_ik=A(:,5); %vector of total expenditures, (NxK)x1
cons_expenditures_ik=A(:,6); %vector of consumption expenditures, (NxK)x1
countryid=A(1:N,1); %vector of country ids, Nx1
va_share_ik=A(:,7+N); %sectoral value added shares of gross output, (NxK)x1
lab_share_ik=A(:,8+N);% labor share of value added, (NxK)x1
cap_share_ik=A(:,9+N); % capital share of value added, (NxK)x1
int_share_iks=A(:,10+N:9+N+K); % intermediate shares, row is using sector and column is input sector, (NxK)xK
%Note: no intermediate shares for ROW (country number 62).

%Parameter vectors
theta_nm=4.32; %non-manufacturing sector TE set to median of manufacturing
h_theta = theta_nm*[ones(2,1);zeros(K_man,1);ones(K-2-K_man,1)]+[zeros(2,1);B(:,2);zeros(K-2-K_man,1)]; %trade elasticities
h_gamma = [zeros(2,1);B(:,1);zeros(K-2-K_man,1)]; %Heterogeneous scale elasticities, non-manufacturing set to zero
rho=1.78; %Upper tier elasticity of substitution


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reshaping data vectors

%Reshape trade matrix to be NxNxK (3 dimensional array)
%rows are exporters, columns are importers, 3rd dimension is sectors
X=zeros(N,N,K);
for k=1:K
    X(:,:,k) = X_ink(1+(k-1)*N:k*N,:);
end
X_ink=X;

%From Gross Output to Value Added
V_ik=GO_ik.*va_share_ik;
V_ik = reshape(V_ik, [N 1 K]);
V_i=sum(V_ik,3);
GO_ik=reshape(GO_ik, [N 1 K]);

%Reshaping share vectors
va_share_ik=reshape(va_share_ik,[N 1 K]);
lab_share_ik=reshape(lab_share_ik,[N 1 K]);

%Intermediate Shares
int_shares=zeros(N,K,K); %First Dimesion is countries, then using sector, then input sector

    for k=1:K
           int_shares(:,k,:)=int_share_iks(1+(k-1)*N:k*N,:);
    end
    int_share_iks=int_shares;

%Reshaping some additional vectors
D_i =reshape(D_i, [N 1 1]);
expenditures_ik = reshape(expenditures_ik,[N,1,K]);
expenditures_nk=permute(expenditures_ik,[2 1 3]);

cons_expenditures_ik = reshape(cons_expenditures_ik,[N,1,K]);
cons_expenditures_nk=permute(cons_expenditures_ik,[2 1 3]);

%Create trade shares
expenditures_nnk=repmat(expenditures_nk,[N 1 1]);
lambda_ink = X_ink./expenditures_nnk; %trade shares

%Create CONSUMPTION expenditure shares
cons_expenditures_n=sum(cons_expenditures_nk,3);
beta_nk=cons_expenditures_nk./repmat(cons_expenditures_n,[1 1 K]);

%Create sectoral domestic trade shares
dom_lambda_ik=zeros(N,K);
for i=1:N
    for k=1:k
        dom_lambda_ik(i,k) = lambda_ink(i,i,k);
    end
end
        
%Create total sectoral exports
exports_ik=zeros(N,K);
for i=1:N
    for k=1:k
        exports_ik(i,k) = GO_ik(i,k)- X_ink(i,i,k);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from Industrial Policy

%Defining output vectors
gfop_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)
gftp_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)
gfip_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)

%Output vectors for Harberger analysis
GO_VA_opt=zeros(K,N);  % Gross output by sector divided by total value added at the optimum
GO_real_opt=zeros(K,N);  %percent change in gross input allocated to sector k when going from optimum to original equilibrium
fs_opt=zeros(K,N);   %foreign expenditures on domestic goods, as a share of value added, at the optimum
fGO_real_opt=zeros(K,N);  % (negative) percent change in effective gross input units supplied to foreigners when going from the optimum to original equilibrium

GO_real_tp=zeros(K,N);  %percent change in gross input allocated to sector k when going from trade policy equilibrium to original equilibrium
fGO_real_tp=zeros(K,N);  % (negative) percent change in effective labor units supplied to foreigners when going from the trade policy equilibrium to original equilibrium

GO_real_ip=zeros(K,N);  %percent change in gross input allocated to sector k when going from industrial policy equilibrium to original equilibrium
fGO_real_ip=zeros(K,N);  % (negative) percent change in effective gross input units supplied to foreigners when going from the industrial policy equilibrium to original equilibrium

GO_VA_initial = zeros(K,N);  % Gross output by sector divided by total value added at the initial equilibrium
fs_initial=zeros(K,N);   %foreign expenditures on domestic goods, as a share of value added, at the initial equilibrium

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from trade and industrial policy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    
%Het TE, Het SE (service 0 SE)
theta=h_theta; 
gamma=h_gamma;
counter=0;

for i=1:N-1 %Excludes ROW
    %Defining data vectors and parameters
    V = V_i(i); %total value added
    X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
    D= D_i(i)/V; %Aggregate trade balance (divided by value added)
    beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
    lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
    V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
    int_share=squeeze(int_share_iks(i,:,:)); %intermedidate goods shares
    va_share=va_share_ik(i,:)';
    %sum(va_share+int_share*ones(K,1))/K
    lab_share=lab_share_ik(i,:)';
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Optimal trade and industrial policy
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=gamma;
    %gamma=zeros(K,1);
    %Solving the model
    model_output=gfop_soe_ces_cap_int_go_solver(V_k,D,X,lambda,beta,int_share,va_share,lab_share,gamma,theta,rho,t_k,s_k,K)
    gfop_hte_0(i,1)=model_output.welfare;
    gfop_error(i,1)=model_output.Z_w;
    
    GO_VA_opt(:,i)=model_output.GO_VA_share;
    GO_real_opt(:,i)= model_output.GO_real; 
    fs_opt(:,i)= model_output.fs_share; 
    fGO_real_opt(:,i)=model_output.fGO_reallocation;
    
    GO_VA_initial(:,i) = V_k./va_share;
    fs_initial(:,i)= X;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Trade policy only
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=zeros(K,1);
    %Solving the model
    model_output=gfop_soe_ces_cap_int_go_solver(V_k,D,X,lambda,beta,int_share,va_share,lab_share,gamma,theta,rho,t_k,s_k,K)
    gftp_hte_0(i,1)=model_output.welfare;
    
    GO_real_tp(:,i)= model_output.GO_real; 
    fGO_real_tp(:,i)=model_output.fGO_reallocation;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Industrial policy only
    %Specifying taxes and subsidies
    t_k=zeros(K,1);
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_cap_int_go_solver(V_k,D,X,lambda,beta,int_share,va_share,lab_share,gamma,theta,rho,t_k,s_k,K)
    gfip_hte_0(i,1)=real(model_output.welfare);
    
    GO_real_ip(:,i)= model_output.GO_real; 
    fGO_real_ip(:,i)=model_output.fGO_reallocation;
    counter=counter+1
end

%Computing gains from optimal trade and industrial policy
gfoip_hte_0 = gfop_hte_0-gftp_hte_0;
gfotp_hte_0 = gfop_hte_0-gfip_hte_0;

%Saving the gains
gains_hte_0=[countryid gfop_hte_0 gftp_hte_0 gfip_hte_0 gfoip_hte_0 gfotp_hte_0];
csvwrite('gains_soe_hte_ces_cap_int_go_0.csv',gains_hte_0);


%Open Economy -FORMULA USING SHARES AT THE Initial Equilibrium
%Create value added shares and foreign sales shares from initial equilibrium.  Create change in sector scale and change in effective labor devoted
%to exports when going from equilibrium with optimal policy to
%laissez-faire equilibruim.  Use these to compute the formula.

%Note: we have a few zero export sectors, for which the formula delivers an NAN
%value for foreign labor reallocation.  Need to replace with zero (which is
%appropriate)

fGO_real_opt(isnan(fGO_real_opt))=0;
fGO_real_tp(isnan(fGO_real_tp))=0;
fGO_real_ip(isnan(fGO_real_ip))=0;

GO_real_opt_tp=zeros(K,N);
for i=1:N
    reallocation_opt =1./(1-GO_real_opt(:,i))-1;   %l_hat -1 when going from initial to optimum
    fl_reallocation_opt =1./(fGO_real_opt(:,i)-1)-1;   % percent change in effective labor units supplied to foreigners when going from the original equilibrium to optimum
    dw_harberger_opt(i,1) =(1/2)*sum(GO_VA_initial(:,i).*reallocation_opt.*gamma+fs_initial(:,i).*fl_reallocation_opt.*-1./(theta+1));% Harberger gains from optimal policy
    dw_harberger_opt_trade(i,1)=(1/2)*sum(fs_initial(:,i).*fl_reallocation_opt.*-1./(theta+1)); %Optimal gains due to trade policy
    dw_harberger_opt_industrial(i,1)=dw_harberger_opt(i,1)-dw_harberger_opt_trade(i,1); %Optimal gains due to industrial policy
    %Now computing Harberger gains from trade policy
    lf_hat_opt = 1./(fGO_real_opt(:,i)-1);  %hat change in foreign labor when going from original equilibrium to optimum
    lf_hat_ip = 1./(fGO_real_ip(:,i)-1); % hat change in foreign labor when going from original equilibrium to industrial policy only equilibrium
    fGO_real_opt_ip=lf_hat_opt-lf_hat_ip;  % change in effective labor units supplied to foreigners when going from the industrial policy only equilibrium to optimum, relative to laissez-faire labor allocation
    fGO_real_opt_ip(isnan(fGO_real_opt_ip))=0;
    dw_harberger_trade(i,1)=(1/2)*sum(fs_initial(:,i).*fGO_real_opt_ip.*-1./(theta+1)); %Harberger gains from optimal trade policy
    %Now computing Harberger gains from industrial policy
    lhat_opt = 1./(1-GO_real_opt(:,i));    %hat change in domestic labor allocation when going from original equilibrium to optimum
    lhat_tp=1./(1-GO_real_tp(:,i));    %hat change in domestic labor allocation when going from original equilibrium to equilibrium with trade policy only
    GO_real_opt_tp(:,i)=lhat_opt-lhat_tp;  %percent change in labor allocation when going from tp only equilibrium to optimum, relative to initial labor allocation
    dw_harberger_industrial(i,1)=(1/2)*sum(GO_VA_initial(:,i).*GO_real_opt_tp(:,i).*gamma);
end
z=[countryid dw_harberger_opt dw_harberger_trade dw_harberger_industrial];
csvwrite('gains_harberger_cap_int_go.csv',z);

%Now saving initial share and reallocation terms to file for comparision
%for case without intermediates
m1 = [countryid GO_VA_initial'];
csvwrite('go_va_initial.csv',m1);
m2 = [countryid GO_real_opt_tp'];
csvwrite('go_real_opt_tp.csv',m2);


